In [ ]:

    

from __future__ import division, print_function

# Third-party
import astropy.units as u
from astropy.constants import G as _G
import h5py
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
%matplotlib inline
pl.style.use("apw-notebook") # comment this out to run at home
importgala.potential as gp
fromgala.units import galactic

import biff
G = _G.decompose(galactic).value


    

In [ ]:

    

def _plummer_density(x, y, z, M, r_s):
    """
    .. warning::

        THIS IS AN INTERNAL FUNCTION -- USE ``plummer_density()`` INSTEAD.

    """
    r2 = x*x + y*y + z*z
    return (3*M / (4*np.pi*r_s**3)) * (1 + r2/r_s**2)**(-5/2.)


    

In [ ]:

    

true_M = 1/G
true_r_s = 1.

x = np.logspace(-2,1,512)
xyz = np.zeros((len(x),3))
xyz[:,0] = x

pot = gp.PlummerPotential(m=true_M, b=true_r_s, units=galactic)
true_pot = pot.value(xyz.T).value
true_dens = pot.density(xyz.T).value
true_grad = pot.gradient(xyz.T).value.T


    

In [ ]:

    

nmax = 16
lmax = 0

Snlm = np.zeros((nmax+1,lmax+1,lmax+1))
Serr = np.zeros((nmax+1,lmax+1,lmax+1))
Tnlm = np.zeros((nmax+1,lmax+1,lmax+1))
Terr = np.zeros((nmax+1,lmax+1,lmax+1))

nlms = []
for n in range(nmax+1):
    for l in range(lmax+1):
        for m in range(l+1):
            nlms.append([n,l,m])
       
for nlm in nlms:
    n,l,m = nlm
    print(n,l,m)
    (S,S_err),(T,T_err) = biff.compute_coeffs(_plummer_density, nlm=nlm, 
                                              M=true_M, r_s=true_r_s, args=(true_M,true_r_s),
                                              epsrel=1E-9)
    Snlm[n,l,m] = S
    Serr[n,l,m] = S_err
    Tnlm[n,l,m] = T
    Terr[n,l,m] = T_err
        
# OR: load from file..
# with h5py.File("/Users/adrian/projects/ophiuchus/data/Anlm_plummer.h5") as f:
#     nmax = f['nlm'].attrs['nmax']
#     lmax = f['nlm'].attrs['lmax']
#     nlm = np.array(f['nlm'])
#     _Anlm = np.array(f['Anlm'])
#     Anlm_err = np.array(f['Anlm_err'])
    
#     ix = np.abs(_Anlm) > np.abs(Anlm_err)
#     nlm = nlm[ix]
#     _Anlm = _Anlm[ix]

# Anlm = np.zeros((nmax+1, lmax+1, lmax+1))
# for (n,l,m),A in zip(nlm, _Anlm):
#     Anlm[n,l,m] = A


    

In [ ]:

    

pl.semilogy(np.array(nlms)[::2,0], np.abs(Snlm.flat/Snlm[0,0,0])[::2])
pl.xlim(0,10)
pl.ylim(1E-6, 1)


    

In [ ]:

    

bfe_dens = biff.density(xyz, Snlm, Tnlm, nmax, lmax, true_M, true_r_s)
bfe_pot = biff.potential(xyz, Snlm, Tnlm, nmax, lmax, G, true_M, true_r_s)
bfe_grad = biff.gradient(xyz, Snlm, Tnlm, nmax, lmax, G, true_M, true_r_s)


    

In [ ]:

    

(true_pot / bfe_pot)[0], (true_dens / bfe_dens)[0], (true_grad / bfe_grad)[0,0]


    

In [ ]:

    

fig,axes = pl.subplots(3, 1, figsize=(8,10), sharex=True)

axes[0].set_ylabel(r'$\rho$')
axes[0].loglog(x, true_dens, marker=None)
axes[0].loglog(x, bfe_dens, marker=None)

axes[1].set_ylabel(r'$\Phi$')
axes[1].semilogx(x, true_pot, marker=None)
axes[1].semilogx(x, bfe_pot, marker=None)

axes[2].set_ylabel(r'${\rm d}\Phi/{\rm d}x$')
axes[2].semilogx(x, true_grad[:,0], marker=None)
axes[2].semilogx(x, bfe_grad[:,0], marker=None)

fig.tight_layout()


    

In [ ]:

    

fig,axes = pl.subplots(3, 1, figsize=(8,10), sharey=True)

axes[0].set_ylabel(r'$\rho$')
axes[0].loglog(x, np.abs((true_dens-bfe_dens)/true_dens), marker=None)

axes[1].set_ylabel(r'$\Phi$')
axes[1].loglog(x, np.abs((true_pot-bfe_pot)/true_pot), marker=None)

axes[2].set_ylabel(r'${\rm d}\Phi/{\rm d}x$')
axes[2].loglog(x, np.abs((true_grad[:,0]-bfe_grad[:,0])/true_grad[:,0]), marker=None)

axes[0].set_ylim(1E-12, 1E0)
fig.tight_layout()


    

In [ ]: